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1 Introduction 

Top quark physics at the LHC is benefiting from huge samples of events that can be used to systemati- 
cally explore many interesting properties of the top quark such as its mass, its production cross-section 
and the structure of its couplings to other Standard Model particles, including the Higgs boson. Many 
of such studies require a good understanding of the kinematics of top quark decay products that are 
affected by QCD radiative corrections. 

The next-to-leading (NLO) QCD corrections to the top quark total decay width and many kine- 
matic distributions are known since long ago [1-3]. The electroweak corrections have been computed 
in Rcfs. [4]. The fully differential computation of corrections to the top quark width at next-to- leading 
order was reported in recent years in Rcfs. [5-7]; it was performed in the context of NLO QCD calcu- 
lations for a number of processes in top quark physics where radiative corrections to the production 
and decay stages are incorporated consistently [5-11]. 

While the NLO QCD description of many hadron collider processes provides quite a reasonable 
approximation, for some processes it may be useful to improve on that and compute the next-to- 
next-to-leading order (NNLO) QCD corrections. Top quark pair production in hadron collisions is 
one of such processes and computations of NNLO QCD corrections to pp — > tt are well underway 
[12-14]. So far these calculations arc performed for stable top quarks and address the total cross- 
section but, eventually, they will mature enough to include also top quark decays in the narrow 
width approximation. This approximation, when done consistently, includes radiative corrections to 
production and decay stages and gives an excellent description of the process pp —> W + W~bb at NLO 
QCD, as was shown explicitly in Ref. [15]. There is no doubt that the narrow width approximation will 
work very well also at next-to-next-to-lcading order, once all relevant corrections to the production 
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and decay stages are accounted for at a fully-differential level. The goal of this paper is to provide a 
fully-differential NNLO QCD calculation for the top quark width. 

The NNLO QCD radiative corrections to the total width of the top quark were computed in a 
number of papers [16-19] and, more recently, NNLO QCD corrections to W-helicity fractions were 
reported in Ref. [20]. However, NNLO QCD corrections to other kinematic distributions in top decays 
were unknown. We note that very recently NNLO QCD radiative corrections to top quark decay at a 
fully differential level were computed in Ref. [21], using a slicing method inspired by soft-collinear effec- 
tive theory [22]. Our computation employs a totally different method and, as such, is complementary 
to Ref. [21]. 

A different phcnomenological motivation for this computation comes from a similarity between 
two processes, t — > b + e + +v and b — > u + e~ + v. This similarity implies that if NNLO QCD radiative 
corrections are known for an arbitrary invariant mass of the lepton pair, i.e. beyond the narrow W- 
width approximation, they can be immediately applied to the calculation of NNLO QCD corrections to 
moments of kinematic distributions in the b — > u semileptonic transition. Since the b — > u decay suffers 
from significant background from b — > c transition, the inclusive rate for b — > u is of little help for 
improving the precision with which the CKM matrix element \V u b\ can be extracted from experimental 
measurements. On the other hand, the knowledge of radiative corrections to moments of kinematic 
distributions with cuts on the electron energy and hadronic invariant mass is useful for constructing a 
theoretical framework to extract \V u b\, see e.g. Refs. [23-26]. The computational framework that we 
develop in the current paper is immediately applicable to the calculation of these moments. 

The other goal of this paper is to contribute to the development of a robust computational scheme 
that will be applicable to any collider process at NNLO QCD. A promising way towards such a general 
scheme was described by Czakon [27, 28] who suggested to combine ideas of sector decomposition [29- 
31] and the Frixione-Kunszt-Signer (FKS) phase-space partitioning [32]. This computational scheme 
was further discussed in Ref. [33] where it was applied to the calculation of NNLO QED corrections 
to Z — > e + e~. In this paper, we employ these methods to compute NNLO QCD radiative corrections 
to top quark decay. From the methodological point of view, top quark decay t — > b + e + + v provides 
an interesting case where the structure of soft and collincar singularities is relatively simple, yet a 
massive colored particle in the initial state leads to a computational complexity 

The paper is organized as follows. In the next Section, we discuss the general set up of the 
computation. In Section 3 we explain how the phase-space for various subprocesses related to top 
quark decays can be parametrized to enable the extraction of infra-red and collinear singularities. In 
Section 4 we discuss the singular limits of various amplitudes required for our computation. In Section 5 
we explain some peculiarities that arise when we try to keep all momenta of external particles in four 
dimensions, yet use dimensional rcgularization consistently. In Section 6 we describe the results of 
the calculation and show some kinematic distributions. In Section 7 we present our conclusions. 
Finally, we give explicit expressions for tree amplitudes employed in our calculation in Appendix A 
and show how to compute one-loop master integrals needed in the constructions of one-loop soft limits 
in Appendix B. 

2 The set up 

As we already mentioned in the Introduction, our goal in this paper is to compute NNLO QCD 
corrections to the fully differential decay rate of the top quark t{pt) — > b(pb) + Ve{p§) + e + {po). 
Throughout the paper, we work in the top quark rest frame. The decay to e + i> is facilitated by an 
emission of a W-boson, which can be either on or off the mass-shell. The interaction of the virtual W 
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boson with the neutrino-positron pair W — > v(p§) + e + (pe) is described by the vertex 

^<5|7l6], (2.1) 

where we use the bra and ket notation for massless spinors with definite helicities. We employ the 
standard notation u + (pi) = = \i), u-(pi) = ) = \i] and note that further discussion of 

spinor-helicity methods can be found in the review Ref. [34]. 

We are interested in the NNLO QCD corrections to semileptonic top quark decay. To compute 
these radiative corrections, we require the following ingredients 

• two-loop virtual corrections to t — > b + v + e + ; 

• one- loop virtual corrections to t — > b + g + v + e + ; 

• tree-level double-real emission contributions to the decay rate, t — > b + g + g + v + e + and 
t^>b + qq + v + e + . 

We note that most of the required contributions are available in the literature. Indeed, two- 
loop QCD virtual corrections to t — > Wb can be obtained from an analytic computation of similar 
corrections to b — > W*u, reported recently by a number of authors [35-38]. One-loop amplitudes for 
the t — > b + g + v + e + process can be extracted from Ref. [6] . Double-real emission amplitudes for 
gluons and quarks are not known in a compact analytic form but can be easily calculated using the 
spinor-helicity formalism. 

However, similar to other NNLO computations, the availability of required amplitudes does not 
make it obvious how to arrive at a physical answer. This happens because of infra-red and collinear 
singularities that need to be extracted from individual contributions before integrating over phase- 
spaces of final-state particles. Until recently, the two methods 1 that were used to extract singularities 
in practical NNLO computations were sector decomposition [29-31] and antenna subtraction [43, 44]. 
Both of these methods were developed further in recent years [45, 46] with the goal to make them 
applicable to more complex processes than previously considered. 

The method that we use in this paper was suggested by Czakon [27, 28] (see also [33]), who 
pointed out that a combination of sector decomposition [29-31] and the FKS phase-space partitioning 
[32] allows to resolve overlapping soft and collinear singularities in a universal way. We note that from 
the technical point of view, the case of top quark decay is particularly simple since no complicated 
phase-space partitioning is, in fact, required. Indeed, singularities in matrix elements occur when 
emitted gluons are either soft or collinear to other massless partons. Since there are only three strongly 
interacting massless particles in double-real emission amplitudes for top decays, there is just one triple- 
collinear limit and no need for angular partitioning arises. The soft singularities are disentangled by 
ordering gluons with respect to their energies and by removing the corresponding 1/2! factor from the 
phase-space. Note that this argument needs to be modified for sub-processes with two quarks in the 
final state and we explain how to do it in what follows. 

We now remind the reader about the general approach to fully differential NNLO QCD compu- 
tations based on sector decomposition and phase-space partitioning [27, 28] (see also [33]). The main 
idea is to express the phase-space for the final state particles through suitable unit-hypercube variables 
that make the singular behavior of scattering amplitudes manifest. Once this is accomplished, it is easy 



1 Other approaches to NNLO calculations are discussed in [39—41]. Some phenomenological applications can be found 
in [42], 
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to construct subtraction terms for real emission processes by using expansion in plus-distributions. To 
be specific, we consider a double-real emission contribution to the differential decay rate of the top 
quark 

rf * } = J dLips \M({p})\ 2 J({p}), (2.2) 

where dLips is a phase-space for t — > b + g + g + e + + v, M({p}) is the decay amplitude and J({p}) is 
a measurement function that depends on momenta of final state particles {p}. We assume that a map 
of the phase-space onto a unit hypercube is found; variables that parametrize the map are denoted as 
{x}. As we will show below, it is possible to re- write Eq.(2.2) as 

Ns W ps J, 

rr , = £ r ( M , 0j pCniM, = jj I _^-F i(e , { x}) (2.3) 

where N$ is the number of sectors, e = (4 — g?)/2 is the parameter of dimensional regularization, N ps 
is the number of independent variables required to parametrize the phase-space, and the function F 
is the product of the scattering amplitude squared, sector-dependent powers of singular variables and 
the phase-space weight Ps^dx}) that may include the measurement function 

F l ({x})=x b ^x b 2 \.x b ^ : \M({p})\ 2 Ps„,({x}). (2.4) 

The important property of these functions F{ is that limits exist at each potentially singular point 

Xj = 0, 

\jmFi(..,x j ,..) = F i (..,0,..) J t:oo. (2.5) 

Xj=0 

The existence of limits Eq.(2.5) allows us to re- write Eq.(2.3) in an explicitly integrable form by 
employing the decomposition in plus-distributions 



r l+ae 



1 r/ N ^ (~ ea )" 
—6{x) + > p- 



n=Q 



ln n (x) 



(2.6) 



Substituting this expansion into Eq.(2.3) and collecting terms of the same order in e, we obtain 
explicitly integrable terms that can be divided into two broad categories, 

1 dxi ^(x i )(F(..,x i ,..)-F(..,0,..)) and j dXi 5(x i )F(..A..). (2.7) 



o o 

It is easy to recognize that the first integral in the above equation contains a "subtraction term" 
and the second integral defines an "integrated subtraction term" , if the language familiar from NLO 
QCD computations is borrowed. It is remarkable that, in the described framework, the subtraction 
terms are obtained automatically, once a proper parametrization of the phase-space is found and the 
expansion in plus-distributions is performed. We also note that any subtraction term depends on the 
function F(x\, ..,xn ) at such values of Xi for which at least one of the singular variables vanishes. 
These kinematic points correspond to universal soft and collinear limits of scattering amplitudes so 
that F(..,Xi = 0,..) can be calculated in a process-independent way. In what follows, we discuss 
essential elements of this construction such as the phase-space parametrization and singular limits of 
amplitudes. 
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Before proceeding to that discussion, we explain some subtleties that arise when one-loop correc- 
tions to the process t — >b + g + e + + is are analyzed within this framework. In this case, the singular 
part of the phase-space is parametrized in terms of two variables, x\ and x%, that describe the energy 
of the gluon g and the relative angle between the gluon and the 6-quark, respectively. The contribution 
to the top quark decay rate is given by the integral 



dxi Ax 2 n 



M 



1 — loop 



M 



(0),* 



(2.8) 



where ellipses stand for other non-singular phase-space factors. While the integrand in Eq.(2.8) sug- 
gests that singularities can be isolated by constructing an expansion in plus-distributions, it is in fact 
not possible to do that right away. The reason is that the one-loop amplitude ■M^Zj^e+j, has branch 
cuts in the limits Xi — > and x 2 — > and we must account for them when constructing the expansion 
of the integrand in plus-distributions. To accomplish that, we write 



[x\x 2 \ 



1 -loop A>(0),* 
Jv h-+bge.+ v Jvl t-+hge+v 



F 1 {x 1 ,x 2 , ■■) + {x\xi) e F 2 (xi,x 2 , •■), 



(2.9) 



where the two functions F\^ 2 are free from branch cuts in x\. 2 — > limits. 

The calculation of real-virtual corrections now proceeds as follows. We insert Eq.(2.9) into 
Eq.(2.8), combine non-integer powers of x\, 2 and expand in plus-distributions. Upon doing so, we 
find a variety of terms that are similar to what we showed in Eq.(2.7). However, in this case the 
integrand is written in terms of Fi(x\,x 2 ) and ^2(^1,^2) and we need a prescription that will allow 
us to compute these functions from scattering amplitudes. To understand how this is done, note that 
Fi, 2 ({x\) may appear in the integrand with one or both of its arguments set to zero, which corresponds 
to soft or collinear kinematics of the final state gluon. To compute F1.2 in those cases, we need to 



investigate 



t-^-bge+v 



(0), 



in soft and collinear limits. Consider first the soft limit, x\ — > 0. 



In this case, the interference of the one-loop and tree matrix elements reads 



lim 2Rc 



M 



1— loop 

t—*bge Jr v 



M 



(0),* 



g 2 s C F S tb (2Re 



t^be+v t^-be+v 



+ g 2 s C A S, 



tb 



M (o) M {o),* 



(2-10) 



where St b and sj: b arc tree and one-loop eikonal factors, respectively. Wc note that in the right hand 
side of Eq.(2.10), the dependence of the one-loop amplitude on the momentum of the soft gluon is 
entirely contained in the eikonal factors. To identify Fi_ 2 at vanishing value of xi, it is important to 
realize that the leading-order eikonal factor Stb is free from branch cuts while the one-loop eikonal factor 

x\ — > 0, the first term on 



Sl b may (and in fact does) contain them. We conclude that, in the soft limit 
the right hand side in Eq.(2.10) matches to Fx(0, x 2 ) and the second term to F 2 (xi — 0, x 2 )(x\x 2 ) 
A similar situation occurs in the collinear limit x 2 — > 0. In this case 



lim 



j^l-loop M (0),* 
JVl t^bge+v JVl t^bge+v 



Pb-Pg 



P, 



(0) 



b^bg 



^1-loop ^(0),* 



p, 



(i) 



b^bg 



t^be+v t^rbe+is 



(2.11) 



where P b ^^ g are the tree- and one-loop splitting functions; these functions can be derived from collinear 
limits of scattering amplitudes described in Ref. [47]. Similar to the soft limit, the tree-level splitting 



function P b ^l bg does not contain any branch cut in x 2 — > limit, while the one-loop splitting function 
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does. Hence, the first term on the right-hand side in Eq.(2.11) is identified with Fi(xi,X2 = 0) and 
the second one with (x\ x<x)~ € F%{xi, X2 = 0). 

If x\ and X2 are both non- vanishing, we can re- write combinations of Fx(x\, X2) and i<2(£i, ^2) that 
appear in the integrand for the rate, through the one-loop matrix element interfered with the tree-level 
amplitude, using Eq.(2.9). The one- loop amplitude for the t — > b + g + e + + v that we require, can be 
extracted from Ref. [6], where Wt production in bg fusion is discussed. 2 For our purpose, the relevant 
amplitudes need to be crossed to describe the t b + g + W + channel. Finally, the evaluation of the 
one-loop amplitude for t — > b+g + e + + v process requires master integrals; they can be computed using 
the program QCDLoops [49]. We found, however, that results of this program become unreliable in 
extreme kinematic limits; for example when the momentum of the gluon in the final state becomes very 
small Eg < 10~ s m t . To overcome this difficulty, we independently implemented all master integrals 
relevant for the computation of the one-loop amplitude for t — > b + g + e + + v using analytic formulas 
provided in Rcfs. [49, 50]. Our implementation allows us to employ quadruple precision and to reach 
very small gluon energies E g ~ lQ~ 12 m t and very small opening angles between the bottom quark 
and the gluon. 

The importance of reaching small energies and angles follows from Eq.(2.7) which shows that 
integration over hypercube variables {x} should, in principle, be performed on the interval x% £ 
[0, 1]. This implies that full matrix elements that contribute to the function F(x) must be evaluated 
at vanishingly small values of x. It is, of course, impossible to compute them at x = and we 
need to introduce additional rcgularization. Following Refs. [27, 28] we provide this regularization 
by discarding such points in {x}-space where the product of singular variables is smaller than a 
user-defined parameter 6. For example, for double-real emission corrections where four variables are 
singular, we require X1X2X3X4, > 5, while for the single-emission contributions we require X1X2 > 5. 
Values of S that we use range from 1CP 7 to 10 -12 , where extremely small values are employed to test 
the independence of final results on the cut-off parameter and the cancellation of infra-red divergences. 
For phenomenological calculations we find that values of S ~ 10~ 7 — 10~ 9 are small enough to provide 
reasonable results for the second order corrections. 

Finally, we note that Eq.(2.7) allows us to compute arbitrary kinematic distributions in a single 
run of the code. We only need to determine the kinematics of the event - which is uniquely defined 
for a given set of {cc}-variables - and assign the relevant weight to a particular bin in a histogram that 
we want to compute. We do this each time the function F({x}) is calculated. Such implementation 
is standard for how subtraction procedure is employed in parton level Monte-Carlo event generators 
both at next-to-leading order (see e.g. Ref. [48] ) and beyond [51]. 

3 Phase-spaces for various sub-processes in top decays 

In this Section, we discuss the phase-space parametrization for the various sub-processes related to top 
quark decays. We start with the leading order phase-space and then proceed to phase-spaces with one 
and two additional gluons in the final state. We also address the parametrization of the phase-space 
for t — > b + qq + e + + v since it differs from the two-gluon case in important details. 

The leading order process that we deal with is t — > b + W*(e + + u), where the VF-boson can be 
off the mass shell. While the W off-shellness is not an important effect for top quark decays, our 
reasons for doing the calculation this way is explained in the Introduction. To account for the W 

2 We use the Fortran77 implementation of this amplitude in the MCFM program [48]. We slightly extend it to allow 
the use of quadruple precision in FortrarrQO. 
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off-shellness in a systematic way, it is convenient to include the Breit-Wigner propagator from the 
scattering amplitude into the definition of the phase-space. We write 

m _ dp 2 w Lipsft ->■ bW*)dUps(W* -» ul+) 

where Lips stands for Lorentz-invariant phase-spaces for indicated final states, p^y is the invariant 
mass of the lepton pair and Tw is the W-dec&y width. We remove the Breit-Wigner distribution by 
changing variables 

Pw — nriw + m w T w taii£, (3.2) 

so that 

J = d C (3 3) 

^ (Pw ~ m w) 2 + m w T w 2irm w T w ' 
The kinematics of the process implies a restriction on the invariant mass of the lepton pair, < p\y < m 2 . 
To account for it, we re- map £ onto the (0, 1) interval by defining £ — £ m ; n + (£ max — £,m\n)%7- We 
obtain 

i(^-4)t^r 2 v = AWd * 7 ' (3 - 4) 



W W I ' '""W w 

where 

jir Cmax £min ^ 777-^ 

Jvbw = -7; 7; — > ?min = - arctan — — , £ max = arctan 



(3.5) 

After generating the invariant mass of the lepton pair using Eq.(3.4) , we obtain lepton and neutrino 
momenta in the W* rest frame in the standard way 

P u,i = ^(i,±m), (3.6) 

where fti = (sin #5 cos 935, sin #5 sin tps, cos #5). Redefining angular variables, 

cos #5 = 1 — 2xs, (p$ = 2-KXg, (3.7) 
we re-write the dilcpton phase-space as 

dLips(J^* -> vl+) = ^£^£ (3.8) 

07T 

We note that we treat the lepton-neutrino phase-space as four-dimensional, even in higher orders 
of QCD perturbation theory where infra-red and collinear divergences require analytic continuation 
of the space-time dimensionality. It is intuitively clear that it should be possible to do so because 
lepton momenta are infra-red safe observables. However, some subtleties arise when the computation 
is organized this way. We discuss them in Section 5. 

Finally, we describe the parametrization of the t —> bW* phase-space. We work in the rest frame 
of the decaying top quark and define the z-axis by aligning it with the momentum of the 6-quark. The 
momenta read 

Pt = m t (l, 0,0,0), Pb = E b (l,n b ), n b = (0,0,1). (3.9) 
The energy of the &-quark is given by 

E b = (3.10) 
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It is easy to find the leading order phase-space for t —> bW* decay. It reads 



Lips(t -+bW) = — [l-Z%) T(e) = 7^T(e), (3.11) 



8ir \ mf J 4Trm 



where i? ma x = E b is the maximal energy of the &-quark for fixed invariant mass of the lcpton pair and 

m = W^ L ^ f%0 =l + 0(e). (3.12) 



T(2-2e 

We can then write the leading order phase-space for t — >• b + e + + v as 

dn be +„ = T(e)7V B w o5T X dx 7 dx 8 dx 9 = ^^^w A _ PvjM dxydxsdxg. (3.13) 
327r^m t 64:iT z \ mf J 

The e-dependent part of the leading order phase-space T(e) can be neglected provided that exactly the 
same quantity is identified and removed from phase-spaces in higher orders of perturbation theory In 
our calculation, we set T(e) = 1 consistently everywhere. 

We now turn to the phase-space which is relevant for the description of the single emission process 
t — > b + g(p-s) + e + + v. The treatment of the W* decay phase-space and the Breit-Wigner factor is 
identical to the leading order case; the new element is the phase-space for the process t — > W* + b + g, 
which we refer to as dLips(£ — > bgW*). This phase-space is split into "regular" and "singular" parts 

dLips(t -> bgW*) = dLips(Q -> bW*) x [dp 3 ], (3.14) 

where Q=p t -p 3 , [dp 3 ] = d^WK 2 ^"^] and 

dLips(Q bW*) = (27T) D 8^(Q-p b -p w )[d Pb }[dp w }. (3.15) 

The momenta that enter the "regular" phase-space read 

p t = m t (l, 0,0,0), p b = E b (l,fi b ), E b = Q2 -Pw (3 16) 

2(Qo -Q-n b ) 



We find 

dLinsfO -> hW*) = ■ 

4^ Q a -Q-n b \E, 
The gluon momentum is parametrized as 

P3 = S3 (1, sin 03 cos (^3, sin 03 sin i^3, cos 03 ), (3.18) 

where gluon energy and emission angles are computed from unit hypcrcube variables 

E 3 = x 1 E max , cos 3 = l-2x 2 , ^3 = 27ra;3. (3.19) 

Expressed through these variables, the singular phase-space becomes 

r, i r(l — e)PS~ e dx'idx2dx3 r , -, .„ „„. 

[*»] - ^F^ 2C W [x?X2] ' (3 0) 

where PS = 16£'^ lax (l — X2) sin 2 <p3. We note that the term [xfa^] in Eq.(3.20) gets combined with 
the matrix element squared \M t ^ b + e ++v+g\ 2 to render the function F(x), that we discussed in the 
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previous Section. It is this factor from the phase-space that makes F(x) finite both in the soft x\ — > 
and in the collinear X2 —> limits. We find it convenient to factor out a term T(l + e)/[87r 2 (47r) _e ] 
per order in perturbation theory and to write the final result for the differential top quark width as a 
series in a s /{2n). We therefore write the phase-space for the process t — > b + g(p 3 ) + e + + v as 



dir 



T( £ )4^4iVBW 



87r 2 (47r)- 



r(l-e) 



r(l + e)r(l - 2e) 



El 



167r2 Qo-Q-n b 



PS" 



, da;i da; 



2 da^da^dirsdzg [a; 2 ^] 



(3.21) 



where Q = p t - p 3 - 

Finally, we discuss the phase-space for the double-real emission process t — > b+g(p 3 ) +g{pi) +e + +i>. 
Similar to the single-gluon case that we just discussed, the phase-space is split into a regular and a 
singular part. The regular phase-space is identical to the single gluon emission case Eq.(3.17), up to a 
change in the definition of the momentum Q = pt — p 3 — P4- The singular phase-space is given by the 
product of single-particle phase-spaces of the two gluons, g{p 3 ) and g(j>i). We choose the momenta to 
be 



p t = m t (l, 0,0,0), p b = E b (l,H b ), p 3 = E 3 (l,n 3 ), p i = E i (l, n i ), 

n 3 = (sin 6 3 cos ip 3 , sin #3 sin tp 3 , cos #3) , H4 = (sin 64 cos ^34 , sin 64 sin (^34 , cos 64 ) 



(3.22) 



with ip 3 4 = (p 3 + tp 4 . We compute the azimuthal angle ip 3 using ip 3 = 2irx 6 . We note that the number 
of particles that participate in the "QCD part" of the process is sufficiently small, so that we can 
embed momenta of top and bottom quarks and of the two gluons into a four-dimensional space-time. 
By momentum conservation, the momentum of the M^-boson is also four-dimensional, but this does 
not imply that momenta of a lepton and a neutrino can be chosen to be four-dimensional as well. We 
will argue in Section 5 that it is in fact possible to do that provided that a spin- correlation part of 
certain splitting functions is modified. 

The correct choice of phase-space variables is crucial for extracting singularities from products 
of matrix elements squared and phase-space factors. For the two-gluon case, we must consider the 
following singular limits: i) one or both gluons are soft; ii) one or both gluons are collinear to a 
6-quark and Hi) two gluons are collinear to each other. The parametrization of phase-space that leads 
to factorization of these singularities was given in Refs. [27, 28] and we now discuss it in the context 
of top quark decay. We note that by ordering two gluons in energy, we remove the 2! symmetry factor 
from the phase-space. We write the singular phase-space as 



dO(<i-3)dO(<i-3) 

[dp 3 }[d P ME 3 2S0 = 25+2e(27r)M _ 2 dy 3 [sin% 3 )] " 

x Kifc] 1- * Mi - %)P Mi - m)V [A(i - A)p 1/2 - £ lm D ?t~ 2e 



(3.23) 



x (2£ max ) 4 - 4£ 0(£i-6)0(z E 



6) d£id£ 2 d77 3 d?74dA. 



Various variables introduced in the above formula parametrize energies and angles of the (potentially) 
unresolved gluons. We use 



E 3 a — E max £i 



1. 



6 



(3.24) 
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and 



?734 



\m - ml 2 



sin 2 tpA = 4A(1 - A) 



\V3 ~ V4\ 2 



D{V3,V4, A) = m + T] 4 - 2 mm + 2(2A - 1)^1^(1 - V3 ){1 - m ). 
The ry-variables are defined as the following scalar products 

2??3,4 = 1 - "3,4 ■ fib, 27734 = 1 - n 3 • n 4 . 



(3.25) 



(3.26) 



Parametrization of the singular phase-space in Eq.(3.23) is still too complicated to extract all singu- 
larities; further decomposition is required. This is achieved by a sequence of variable changes that we 
describe below following Rcfs. [27, 28]. Specifically, we split the phase-space into a sum of five different 
sectors. In each of these sectors, the relation between "old" variables and "new" variables {x} 

reads 

Z3Z4 

V4 = x 3 , (3.27) 

m = x 3 , 

2J' V4 = X3 - 

We also write A = sin 2 {ttx^/2). This change of variables introduces a factor of 7r in the normalization 
of the phase-space that is included in the expressions below. The integration region for x§ is x§ G [0,1]. 
The full phase-space for t — > b + g + g + e + + v is given by the sum of five (sector) phase-spaces 



sector 1 : 




= a>i, 


6 


= xix ma _ x x 2l m = 


X3, V4 


sector 2 : 


a 


= Xi, 


6 


= XiX max X 2 , T} 3 = 


xs, 774 


sector 3 : 


& 


= Xx, 


6 


= X\X max X 2 x^ rj 3 


X3X4 

2 


sector 4 : 


6 


= SCl, 


6 


= xix max x 2 , r/ 3 = 


X 3 X±X 2 

2 


sector 5 : 


6 


= £1, 


6 


— X±X rnax X 2 ^ Tj 3 — 


x 3 fl- 



dn 



(3.28) 



Each of the five phase-spaces is written as 



dn iLe + , = n °™ >< ps»,iPsr e n dx * x n 



da;. 



5=1 Ij 



b< 4) &<*> 6« 6« 



(3.29) 



Below we present the functions PS lUi i, PSi and the exponents 0^=1 4 and 6!-^ 4 for each of the 



sectors. We note that the normalization factor is common to all sectors; it reads 

2 



Norm = T(e) 
We also note that we can write 

1 



rq + e) 

(47r)- e 87T 2 



PS. 



8 ^ 2 Qo~Q-fi b 



PS, 



New 



PS, = 1024£ b 2 i; i 2 iax A(l - A)< ax (l - x 3 ) sin" (^3) PS*. (3.31) 



r(l-e) 



r(l + e)r(l - 2e) 



(3.30) 
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The expressions for the exponents and the phase-space factors for each of the five sectors read 3 



sector 1 : {a\ = 4, 02 = 2, a 3 = 2, a 4 = 1}, {b\ = 4, b 2 = 2, b 3 = 2, 64 = 1}, 



2 



2iV 1 (x 3 ,f ,A)' 2iV 2 (x 3 ,f ,A) ' 

sector 2: {ai = 4, a 2 = 2, a 3 = 2, a 4 = 2}, {61 = 4, b 2 = 2, fc 3 = 2, fe 4 = 2}, 

4^(^,1- f,A)' 47V 1 2 (.x 3 ,l-f ,A) ' 

sector 3 : {ai = 4, a 2 = 2, a 3 = 2, a 4 = 3}, {61 = 4, fe 2 = 2, fe 3 = 2, fe 4 = 3}, 

(1-^) (1-f) 2 ( 3 - 32 ) 
2^3^, A)' rS 2iV 1 2 ( a; 3 ! f ,A) ' 
sector 4 : {a 4 = 4, 02 = 3, a 3 = 2, a 4 = 1}, {61 = 4, b 2 — 3, & 3 = 2, £> 4 = 1}, 

pc 2 / pc _ y± 2 1 \r 2 ) 

* w 2N 1 (x 3 ,^,\y ° 2JV?(x 3 , 2 f a ,A) ' 

sector 5 : {a 4 = 4, a 2 = 2, a 3 = 2, a 4 = 2}, {6 4 = 4, fe 2 = 2, fe 3 = 2, fe 4 = 2}, 

ps 1 po_ a-f)(i-^(i-f)) 

^ 4^(^,1- f,A)' 47V 2 (x 3 ,l-f ,A) ' 



The function N\ is given by 



7V 1 (x 3 ,x 4 , A) = 1 + 14(1 - 2a; 3 ) - 2(1 - 2X)y/x 4 (l ~ x 3 )(l - x 3 x 4 ). (3.33) 

It is straightforward to use momentum and phase-space parametrizations described above to extract 
singularities from matrix elements squared in a factorized form, thereby facilitating the integration 
over the full phase-space with the help of plus-distribution expansions. We will briefly discuss this in 
the next Section. 

We note that the phase-space parametrization that we just described resolves all the singularities of 
the matrix elements for the final state with a &-quark and two gluons. The process with a 6-quark and 
two additional masslcss quarks in the final state has a simpler singularity structure than the two-gluon 
case, so our phase-space parametrization should be applicable to the quark case as well. However, 
the quark case is different in that quark amplitudes are, in general, not symmetric with respect to 
q and q permutation and also they are less singular. It is therefore useful to consider bqq final state 
separately and design a parametrization of the phase-space for this case. We start by listing kinematic 
configurations that lead to non-integrable singularities in the amplitude t — > b + q(p 3 ) + g(p 4 ) + e + + v. 
From the point of view of the singularity structure it is sufficient to understand the singlet contribution 
where the flavor of the 55-pair is different from b. The case of bbb final state then easily follows. The 
singular configurations that we need to consider are i) both q and q are soft; ii) both q and q are 
collinear to the 6— quark and Hi) q and q are collinear to each other. Analyzing these singular limits, 
we conclude that the number of sectors that we need to consider for bqq final state is four. Momenta 
are parametrized as in the gluon case. The first two sectors are 

• sector 1: E 3 = XiE max , E4 = XiX 2 x max E max , rj 3 = x 3 , r\<± = x 3 (l - x 4 ), 

• sector 2: E 3 = x\E max , E 4 = XiX 2 x max E max , r] 3 = x 3 (l - x 4 ), r\<± = x 3 , 



3 We suppress the sector label everywhere in the equations below. 
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and the other two are obtained from sectors 1 and 2 by permuting p 3 and p±. Finally, all quark sectors 
share the same phase-space. It reads 



r(l + e) 1 



(47r)- e 8tt 2 

1 -£'6-E'max a; max x ^ 



2 



AW 



r(i 



8tt 2 (Q -Q 3 ) JVi(aj 3 ,l-X4,A) 
4S7 4S74S7dx 2 d, 5 d, 6 dx 7 dx 8 dx 9 [,fx^ 2 ], 

Xi 



r(l + e)r(l-2e), 
PS~ e (3.34) 



where 



ps = L-Bg -" A (l - A)(l - x 3 )(l - x 4 )(l - x 3 (l - x 4 )) sin 2 y 3 . (3.35) 

N{(x 3 , 1 - £4, A) 

It follows from the phase-space parametrization that we do not intend to construct the plus-distribution 
expansion with respect to the variable x 2 . This is done intentionally since x 2 — > corresponds to a 
single-soft limit which is not singular for the qqb final state. 

4 Singular limits 

The phase-space parametrization that we discussed in the previous Section is optimized for extracting 
the singular behavior of scattering amplitudes. This is done by exploiting universal soft and collinear 
limits of these amplitudes. The goal of this Section is to describe the limits that we require for the top 
quark decay. We begin with the discussion of the relevant soft limits and then consider the collinear 
ones. 

Soft singularities can only occur if gluons become soft. 4 In general, emissions of soft particles lead 
to color correlations but because in the top quark decay the number of colored particles is small, all 
color correlations simplify and the relevant soft limits can be written in a compact form. For example, 
the single soft limit of the t — > b + g(p 3 ) + e + + v amplitude squared reads 

\M(t, b, g)\ 2 w g 2 s C F \M(t, b)\ 2 S tb { P3 ), (4.1) 

where the eikonal factor Stb is 

2p t ■ Pb 

{Pt ■ P3){pb ■ p 3 ) (pt-p 3 y 
Note that for a single- gluon emission process t — > b + g(p 3 ) + e + + v the function F{x\ ,x%) is given by 

F{x 1 ,x 2 )=x 2 1 x 2 \M{t 7 b,g)\ 2 7 (4.3) 

where x\ and x 2 parametrize the gluon energy E 3 = E mayL x\ and the relative angle with respect to 
the 6-momentum, Pb ■ p 3 = '2EbE 3 x 2 - We have seen that we need F(x\ — 0, x 2 ) for the computation 
and we now show how to obtain it from the soft limit Eq.(4.1). The key point is that the dependence 
on the gluon momentum in the right hand side of Eq.(4.1) is entirely contained in the eikonal factor. 
Hence, to find F(x\ = 0,22), we need 

lim xlx 2 Stb(p 3 ) = E m l x (l - x 2 ). (4.4) 

x 1 — >0 



S tb{P3) = jz Z \7Z — TT ~ 71 — ( 4 - 2 ) 



4 In the process t — >b + q + q + e + + v there is also a soft singularity when energies of both q and q vanish. We do 
not discuss this case as it follows from the soft limit of t — > b + g + e + + v. 
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While the single soft limit provides a very simple example, we note that the same approach is used 
to compute limiting values of integrand functions F also for more complex cases, such as limits of 
double-real emission amplitudes squared and limits of one-loop contributions to single-real emission 
processes. We will not discuss these derivations; instead, we limit ourselves to the presentation of the 
relevant ingredients in what follows. 

The double soft limit of the amplitude squared for t — > b + g(p3) + g(pi) + e + + v is only slightly 
more complicated than the single soft limit. We can write it as 

\M(t,b,g 3 ,g 4 )\ 2 ttg*\M(t,b)\ 2 C F (C A (2S tb (p 3 ,p 4 ) ~ S tt (p3,Pi) - S bb (p 3 ,pi)) 

(4.5) 

+ C F S t b(p 3 )S t b(pi)) ■ 

The double-emission cikonal factors depend on the masses of the emitters; they can be extracted from 
Refs.[28, 52]. For completeness, we present them below for the specific cases that we need to describe 
top quark decays. In the massless case p\ = 0, we find 

a , \ nfi s (Pb ■ Ps) (Pb ■ Pi) fA - 

dbb(P3,P±) = 2(1 - e)- 7,, (4.6) 

[P3 ■piYiPb-pzi) 1 

where 7534 = p 3 + Pi- In the massive case p\ = m 2 , the double eikonal factor reads 

(pt ■ Pz) (pt ■ Pi) m\ 



Stt(P3,Pi) =2(l-e)- , 

(P3 ■ Pi) (Pt ■ P3i) (Pt ■ P3) (Pt ■ Pi) (Pt ■ P3i) 

3 mf m 2 I ' p t ■ p 3 p t ■ pi 



+ -- 7-- — -2 I + 



(4.7) 



2 (P3 ■ Pi) (Pt ■ P3) (Pt ■ Pi) (p 3 ■ pi) (p t ■ P3i) \ Pt ' Pi Pt ' P3 

Finally, we give an expression for the eikonal factor Stb(p3,Pi), where p 2 = m 2 and pi = 0. It reads 

a , s /, s (Pb -Pi) (Pt -P3) + (p b -P3) (Pt -Pi) , Pb-Pt 
i>tb(P3,Pi) = (1 - ej- 



(P3 ■ Pi) (Pb ■ P3i) (Pt ■ P3i) p 3- p i V (Pb ■ P3) (Pt ■ Pi) 

1 1 / 3 + (Pb-Pi) (Pt-Ps) + (Pb ■ Ps) (Pt ■ Pi) 

(Pb ■ Pi) (Pt ■ P3) (Pb ■ P3i) (Pt ■ P3i) V 2 (Pb -P3) (Pt -Pi) 2(Pb -Pi) (Pt -P3), 
(Pb ■ Ptf (Pb ■ P3) (Pt ■ Pi) + (Pb ■ Pi) (Pt ■ P3) 



(Pb ■ Ps) (Pb ■ Pi) (Pt ■ P3) (Pt ■ Pi) V 2 (Pb ■ P3i) (Pt ■ P3i) 

(Pb-Pt) Pb-P3i 1 



2 (Pb ■ P3) (Pb ■ Pi) (pt ■ P3) (Pt ■ Pi) Pt ■ P3i 4 (p 3 ■ pi) (p t ■ p 3 ) (pt ■ Pi) 
1 / (Pb-P3) 2 (Pb-Pif \\ 



2 (P3 ■ Pi) (Pb ■ P3i) (Pt ■ P3i) \ (Pb ■ Pi) (pt ■ Pa) (Pb ■ Pa) (pt - Pi) J ) ' 

It is straightforward to check, using momenta paramctrization in terms of x- variables for various types 
of sectors for double-real emission described in Section 3, that all overlapping singularities get resolved 
in the sum of the eikonal factors Eq.(4.5) and that the subtraction terms for the soft limits of the 
double-real emission process can be calculated. 

In order to subtract soft singularities from one-loop amplitudes for t — > b + .9(^3) + e + + u, we need 
the expression for the eikonal current to one-loop. The massless result for the one-loop soft current 
was obtained in Ref. [53]. The one- loop massive soft current was computed in Ref. [54]. We have 
re-calculated this soft current for the top quark decay into a bottom quark and present the result of 
the calculation below. This result is identical to Ref. [54] but it is written in the kinematic region 
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which is relevant for the top quark decay, so that no analytic continuation is required. The soft limit 
of the interference of one-loop and tree-level amplitudes is written in Eq.(2.10). The one-loop eikonal 
factor reads 

3^ (pa) = g 2 C F S tb (p 3 )2Re [ltb(pt,Pb,p 3 )] , (4.9) 

where 



{.Pt ■ Pb) 

™%(pb -Pz) 2 ~ 2(Pt -Pb)(Pt -P3){Pb -Pz) 



J tb = ~ -271 ~ ^ 7uZ ~ UZ Z7T7Z ^^ { 2 (;Pt ' P^fiPb ' Ps)h 



+ 2(p t -p 3 )(Pt 'Pb)(p& '^3)^2 + 2(p h -_p 3 ) [ip t -Pb){Pt -Pz) - m 2 {Pb -P3)] h (4-10) 

■P3)(Pb ■ 

(pt ■ Pb) 



A\t \ I \( \ 2/ s21 (Pt ■ P3)(Pb ■ P3) T 

4 [(pt ■ Pb)(pt ■ P3)(Pb ■ P3) - m t {pb ■ P3) \ ] ; h 



The Feynman integrals Ii 4 in Eq.(4.10) are defined as 

[ d d k 1 r d d k 1 

1 ~J i{2n) d {k+p 3 ) 2 k 2 {2p t -kY 2 ~ 7 z(2tt)^ fc2( 2pt . k)(2 Pb -k + 2p b - p 3 ) ' 

/ d d fc 1 /• d d fc 1 

h = / T77— 7^7^77- — 7-777- — r — ^- — -ti J 4 



i{2TT) d k 2 {2p b -k){2pt-k-2pfp 3 y J i(2Tr) d (k-p 3 ) 2 k 2 (2p b -k)(2pfk-2pf P3 )' 

We note that each propagator implicitly includes -MO-prescription that unambiguously specifies the 
relevant analytic continuation of the result. The integrals are computed to be (see Appendix B for 
details) 

r(-e)r(2e)m 2e T(l - 2e)T(e)T(2e)m^ 2e e 2me T(-e)T(2e)m 2e 



h = — — — - — - — — Io = — - — U 



I 



4 



(47r)^/2( St3 )l+2 £ ' ( 47r )d/2( S63 )2e( St6 )l-2e ' - (4^/2(^)2^- 

r(-e)r(l + 2e)mf \ „ nie r(l + e)r(-2e) 



(47r)^/2( a;2 ) e ( St 3) 1 + 2e S6 3 
1 

r 



r(i- c ) 

F 2 i(-e, -e, 1 - e; 1 - x 2 ) 



-F 21 (- e) l + e,l-e;l-x 2 ) (4.12) 



where Sy = 2^ -pj and x' 2 = (f7i 2 Sfj3)/(stbSt3)- We emphasize that products of all four-momenta are 
positive-definite since we write the results for the soft current in the correct kinematic region. We note 
that to use the results Eq.(4.12) in our plus-distribution expansion, we should extract all the relevant 
branch-cuts from the hypergeometric functions. 

Wc now turn to the discussion of collinear singularities. It is well-known that in the collinear limit 
any matrix element squared factorizes into a splitting function that is, in general, spin-dependent, and 
a reduced squared matrix element that depends only on the sum of momenta of collinear particles. For 
the top quark decay, we need to describe such collinear splitting processes as b b + g, b — > b + g + g, 
g —t g + g, b — >• b + q + q, b + bb and g ~ > qq. The corresponding splitting functions can be found in 
Ref. [52]. Many splitting functions that we need for top quark decay do not contain spin-correlations. 
The corresponding factorization formula reads 

\M\ 2 ~ s ^P({z})\M\ 2 , (4.13) 

where M is the reduced amplitude, s c is the scalar product of collinear momenta that vanishes in the 
exact collinear limit and P({z}) is a splitting function that depends on the set of energy fractions {z}. 
This is an easy case to deal with since, in this situation, the collinear limit depends on the reduced 
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amplitude squared which can be computed in the standard way. However, there are two splitting 
functions that we have to deal with (g — ¥ gg, g — > qq), that do contain spin correlations. In that case, 
the collincar limit is more complicated 



\M\ 2 ~ s- 1 P" v ({z})M"M v '*. 
Any spin-dependent splitting function P A " y (z) can be written as 

P(z) = -g^P 1 (z,e) + n^P 2 (z,e), 



(4.14) 



(4.15) 



where is a four-vector that specifies how the collincar direction is approached. To compute the 
right hand side in Eq.(4.14), it is useful to re- write it through helicity amplitudes. This is easy to 
do since gauge- invariance of scattering amplitudes allows us to replace g^ v and n^n" in Eq.(4.15) by 
rank-two tensors constructed from polarization vectors 



9^ 



e„(A)e:(A), n^n" 



E 

A 1= ±,A 2 =± 



n-e(X 1 )n-e(X 2 )} e*(AiK(A 2 ). 



(4.16) 



We note that, since we work in dimensional regularization, it is important to extend the sum over 
four-dimensional gluon hclicities in Eq.(4.16), to include e-dimensional polarization vectors as well. 
It is possible to do so by introducing "scalar" polarizations of gluons that point along Cartesian 
components of the e-dimcnsional sub-space of the full <i-dimensional Minkowski space. In Section 5 we 
will return to the issue of spin-correlations, to discuss a consistent implementation of angular variables 
in our phase-space parametrization. 

While the majority of collincar limits is completely determined by the known splitting func- 
tions [52], there are a few cases when further computations are required. This happens for exam- 
ple in the strongly-ordered limit of the b — > g{p 3 )g{p4)b splitting function. The strongly-ordered 
limit occurs when the two gluons are much more collinear to each other than to the 6-quark, i.e. 
Pa • Pi "C Pb • Pa + Pb • Pa- In this case the full matrix element factorizes as 



\M(t,b,9 3 ,gi)\ 2 « To U9 9s ^ 9 ,P s 9 ; g JM(t,b)\ 2 . 

(2p 3 ■ p A ) (2p b ■ p 3 + 2p b ■ pi ) J3yi 



(4.17) 



We derive the following result for the strongly-ordered limit of the b — ¥ g{pz)g{pA)b splitting function 
"4Az 3 (l-z 3 )z& , 1 + (23(1 - z 3 )(l - z b )) 2 + z 2 b - 2z 3 (l - z 3 )(l + z 2 b )' 



P: 



-8C'fCa 



8C fG a 



1 



(1 



- Zb 
Zb ) 



(1 - z b )z 3 (l - za) 



23(1 - z 3 ) 



2:3(1 - za) 



4Az 3 (l - z 3 )z b 
1 - z b 



(4.18) 



e + 0(e 3 



where z 3 = E 3 / '(E 3 + E4), z b = E b / (E 3 + E A + E b ) and A is one of the phase-space variables introduced 
in Section 3. 



5 How to treat spin-correlations consistently in the NNLO computation 

As we pointed out earlier, it is plausible, but not entirely obvious, that one can consistently treat 
positron and neutrino momenta as four-dimensional in higher-order perturbative computations. The 
goal of this Section is to elaborate on this point. We will discuss it in a slightly simplified setting, by 
considering top quark decay into a polarized VK-boson, ignoring VK-decays into leptons. 
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We want to treat the polarization vector of the W -boson as a four-dimensional vector. In addition, 
there are other vectors that we would like to treat as four-dimensional; they include the W-momentum, 
the top quark momentum and the 6-quark momentum. Assuming that we are in the top quark rest 
frame and making use of momentum conservation, we conclude that we can embed the momentum of 
one additional gluon to the four-dimensional space, but that the second gluon must have additional e- 
dimcnsional components. This means that the phase-space parametrization for single-gluon emission 
processes that we discussed in Section 3 is perfectly valid but the phase-space parametrization for 
double-real emission computations must be extended. To this end, we parametrize momenta and the 
VF-polarization vector as 

p t = (m t) 0;0), p b = E 6 (l,5,l;0), p w = (E w ,p x wl 0,p y w ;0), e w = (e&\ e$ , e^, e$; 0), 

(5.1) 

p g3 = E g3 (1, sin # 3 cos ip 3 , s'm8 3 sin ip 3 cos a, cos #3; sin 9 3 sin <p 3 sin a) , 

and determine p gi using momentum conservation. A component of a four- vector in Eq.(5.1) shown 
after the semi-column is the e-dimensional component. We now perform a rotation in the y — e plane, 
to remove the e-dimensional part of the gluon momentum p g3 . From momentum parametrization of 
other particles it follows that this rotation does not induce e-dimensional components for any of them. 
On the other hand, such e-components appear in the polarization vector ew of the W-boson after the 
rotation. 

To understand in which cases we require the e-dimensional components in Fy-polarization, we 
consider singular limits. These limits can be divided into two groups. Some limits lead to expressions 
where squares of reduced matrix elements appear. The reduced matrix elements squared can be either 
leading order or next-to-leading order. In both cases, because of the choice of the momenta above, we 
conclude that 

\M\l ^ {p t ,Pb,pw^w{a)) = \M\l o ^ hO (pt,Pb,pw,£w(0)), (5.2) 

where a parametrizes the e-dimensional angle and a = corresponds to the four-dimensional limit, 
see Eq.(5.1). We note that Eq.(5.2) is valid because through next-to-leading order none of the vectors 
contains y-components that can lead to an a-dependence through products with ew(a). Therefore, 
in all limits that are proportional to the reduced matrix elements squared, we can neglect the depen- 
dence of the W-polarization vector on the extra-dimensional angle a, provided that the normalization 
condition e^ = — 1 is maintained. 

A different situation occurs when spin correlations appear in the singular limits. In the top quark 
decay case, there are no triple-collinear limits which have this feature and we only need to focus on 
doublc-collinear limits. As an example, we take the case when gluon g 3 and gluon are collinear. 
The collinear limit is written as 

\M t ^ bWg3gi \ 2 ~ P^MW = (-.9 M „Pi(z,e) + P 2 (z,e)n^)M^ bWg3i MU bWg ^ (5-3) 

and n M is the vector that describes how vectors p g3>gi approach the collinear direction. The matrix 
element M^ bWg contains the polarization of the W^-boson and depends on momenta of t, b and 
W. Since none of these momenta have y-components, we can remove the a-dependence from the 
matrix element by rotating all vectors in the y — e plane. Since, in general, the vector does have 
a j/-componcnt, such a rotation induces an a-dependence in the splitting function. The vector 
becomes 



-> ±v r \n fJ - + Vl - A (cos an M + sin an£ ) , 



(5.4) 




E, [GeV] E, [GeV] 

(2) 

Figure 1. Distribution of the second order coefficient dF t in the positron energy. Left pane: double-real and 
real-virtual contributions. Original histograms are shown with solid lines; Legendre polynomial fit is shown 

(2) 

with dashed lines. Right pane: distribution of dF t in the lepton energy. See text for details. 



where fz M = [cos #3 cos 933, cos $3 sin tp 3 , — sin# 3 ; 0] , = [— sin ip 3 , cos ips , 0; 0] and n£ = [0, 0, 0; 1]. The 
two possible signs in front of the first vector are for different sectors. Since we do not need the a-anglc 
for any other limit, it should decouple from the point-by-point subtraction. This allows us to average 
over a and modify the spin-correlated part of the splitting function by including "integrated" terms. 
Using the two non-vanishing averages 

(cos 2 a) = ^ 2a " ) = T~2e ^'^ 

we find the replacement rule 

n^n u -> XfPn" + ^ ~ ^ {n»h v - n£n%) . (5.6) 

Applying it to Eq.(5.3) and treating all the vectors in the matrix element Ad, including W-polarization, 
as four-dimensional, we get an easy way to consistently combine dimensional regularization with four- 
dimensional handling of all the relevant vectors that we have in the process. Finally, we note that it 
is obvious how to generalize this discussion to account for the W^-boson decays to a lepton pair. 



6 Results 

In this Section we show some results of our computation. We note that our goal in this paper is not 
to provide extensive phenomcnological studies of top quark decays but to illustrate the capabilities of 
the computational method. The results that we present below are selected accordingly. 

Wc remind the reader that wc consider the top quark decay into hadronic final states and a single 
lepton pair t — > b + Xhadr + e + + v. We use the on-shell renormalization scheme for the top quark and 
the MS renormalization scheme for the strong coupling constant in a theory with five massless flavors. 
We set the renormalization scale to the value of the top quark mass fi = m t . Results for \i ^ m t can 
be immediately obtained by using renormalization group equations, so we do not present them here. 
In our computation, we account exactly for the off-shellness of the VF-boson. The reasons for doing 
that are explained in the Introduction. 

To show the quality of our numerical computation, we isolate all the different color structures that 
contribute to the final result 



1 + Ci^Aj + (g) 2 (C 2 F A A + C F C A A NA + CpTHmAi + C F T R A h 



0{al). (6.1) 
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The different contributions to second order corrections were computed in a number of papers a decade 
ago for top quark decays to an on-shcll W boson. Since ^A,NA,i,hi as denned in Eq.(6.1), are relative 
corrections, we expect that the importance of the VF-off-shcllncss effect is reduced. We use the analytic 
computation of Refs. [18, 19] performed as an expansion in m\y /m t , to infer numerical values for the 
coefficients. We employ the input parameters m t = 172.85 GeV and mw — 80.419 GcV to obtain 
[18, 19] 

C F T R Ai = 7.978, C F T R A h = -0.1166, CpA A = 23.939, C F C A A NA = -134.24. (6.2) 
We compare these results of the analytic computation [18, 19] with our numerical results. We find 5 
C F T R Af um = 7.970(6), C F T R A^ m = -0.11662(1), 

C F A^ m = 24.38(25), C F C A A n ^ = -133.6(4), ^ 

where integration errors are explicitly shown. By comparing Eq.(6.2) and Eq.(6.3), we find good 
agreement between results of analytic and numerical computations. The level of agreement is rather 
striking since there are large cancellations between various contributions to the NNLO final result. To 
illustrate this point, we show individual contributions to the abelian coefficient CpA A um with their 
integration errors. They read 

C F A A>Kn = -127.12(23), C|A A , RV = 60.24(6), C 2 F A AyY = 83.0154(5), 
Cf A A , ren = 16.513(12), C 2 F A A . m = -8.227(2), 

C F A A = C f Aa,rr + C F Aa,rv + C^Aa.vv + C F A A ^ cn + C F A Ajq q = 24.38(25), 

where we include double- real (RR), real- virtual (RV), two-loop virtual (VV), renormalization and 
quark-pair corrections. Comparing individual contributions with the final result for A^, we observe 
very significant cancellations. It is therefore quite remarkable that we can maintain 0(1%) error on 
the NNLO coefficient. 

We now turn to the discussion of kinematic distributions. To facilitate this discussion, we write 
the differential decay rate as 

dr t = drf + (£) dif > + (g) 2 dlf > + 0(a% (6.4) 

(2^ (2) 

and plot dY t jdx in Figs. 1, 2 for various observables x. In particular, in Fig. 1 we show dT t /dEi, 
where E\ is the positron energy. The final result for NNLO QCD corrections arises as a result of 
large cancellations between various contributions, most notably the double-real and the real-virtual 
ones. To illustrate this point, we show these contributions individually in the left pane and the 

(2) 

sum of all contributions to dT t /dEi in the right pane. We note that at smaller values of E\ the 
final result is about ten times smaller than double-real and real-virtual contributions separately. At 
this level of cancellations, bin-bin fluctuations that are barely seen for double-real and real-virtual 
contributions individually (left pane in Fig. 1) become a serious issue. To overcome it, we choose to fit 
individual contributions to the sum of Legendre polynomials and then add fitted results rather than 

(2) 

original histograms. Upon doing that, we find that our final result for dT t /dEi exhibits a re-assuring 
pattern. Suppose, we account for Nl Legendre polynomials in the fit. If Nl is too small, we find 

(2) 

that the resulting dF t /dEi is not stable against including more polynomials into the fit. If Nl is too 

5 All numerical results presented in this Section are obtained using the adaptive Monte-Carlo algorithm VEGAS [55] 
as implemented in the CUBA library [56]. 
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Figure 2. Left pane: distribution of the second order coefficient dT t in invariant mass of the positron 

(2) 

and the hardest jet. Right pane: distribution of the second order coefficient dl^ in the opening angle of the 
positron with respect to the VK-direction of motion, in the VK-rest frame. See text for details. 



large, we fit bin-bin fluctuations and do not gain anything. However, we find that there is a range of 
intermediate values of Nl that we can use in the fit so that, on one hand, our final result for dT^' /dE t 
does no depend on the exact value of Nl and, on the other hand, the resulting distribution is smooth. 
Distributions shown in the right pane of Fig. 1 and in Fig. 2 are obtained following this procedure. 

In the left pane of Fig. 2 we show NNLO QCD contributions to the kinematic distribution in 
the invariant mass of the positron and the hardest (in energy) jet in the event. The jet here is 
defined with the lepton collider fc^-algorithm where the distance between two partons i and j is given 
by Uij = mm (^ 2 / m t j /m 2 )(1 — cosdij). The relative angle Oij is defined in the top quark rest 
frame. For numerical computations, we take yij =0.1. In the right pane of Fig. 2 we show NNLO 
QCD correction to the kinematic distribution of the positron polar angle defined in the VF-boson rest 
frame, relative to the direction of motion of the VF-boson 6 . This distribution is interesting because it 
allows us to determine helicity fractions of the VF-bosons in top decays. Indeed, to all orders in QCD 
perturbation theory, the decay rate can be written as 

= - A sin 2 0,r L + 1 (i + cos etf r + + \ (i - cos erf r_ . (6.5) 

d cos 0i 4 8 8 

The widths Ti,r± define partial decay rates into polarized M^-bosons. The helicity fractions are con- 
structed from partial widths as F± ! l = T±x/^t, where T t = F + +T_ + Tl. Our result for dr t /dcos#z 
shown in Fig. 2 allows us to compute the NNLO QCD corrections to the helicity fractions. Upon doing 
so, we find good agreement with similar results presented in Ref. [20]. For example, by fitting the 
angular distribution shown in the right pane of Fig. 2 we find the NNLO QCD contributions to helicity 
fractions 7 [5F L , 5F- , 5F+] = [-0.0022(1), 0.0021(1), 0.0001(1)]. These numbers should be compared to 
the results of analytic computations reported in Ref. [20], [5F L ,6F-,SF+] = [-0.0023,0.0021,0.0002]. 
A good agreement between the two results is obvious. 

7 Conclusions 

In this paper we described a computation of NNLO QCD corrections to semileptonic decays of the 
top quark at a fully-differential level. We have used a framework described in Rcfs. [27, 28, 33] that 

6 The momentum of the VF-boson can be determined from the momentum of the recoiling hadronic system in top 
decay. 

7 The exact definition of the helicity fractions and values of a a used to obtain these results can be found in Ref. [20]. 
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combines sector-decomposition and Frixione-Kunszt-Signer phase-space partitioning to disentangle 
overlapping kinematic regions that may lead to soft and collinear singularities. We have shown that 
this technique works in the same way as parton level Monte Carlo integrators at next-to-leading order, 
so that weighted events can be generated and a variety of kinematic distributions can be computed in 
a single run of the program. 

There is a number of phcnomenological applications of the calculation reported in this paper that 
we envision. Eventually, it will be combined with the perturbative QCD description of the top quark 
production process, to provide NNLO QCD results for pp —> ti — >• W + W~bb in the narrow width 
approximation. It can also be used to investigate the dependence of various observables studied in 
top quark decays, for example the helicity fractions, on kinematic cuts applied to top quark decay 
products. A minor modification of our calculation turns it into a NNLO QCD result for the most 
important contribution to the single top production process in hadron collisions. Finally, as explained 
in the Introduction, our calculation can be adapted to explore the impact of NNLO QCD corrections to 
inclusive charmless semileptonic decays of &-quarks on the determination of the CKM matrix element 
\Vub\- We plan to return to these applications in the near future. 
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A Tree-level amplitudes 

In this Appendix we collect formulas for helicity amplitudes for top quark decay processes t — > b + 
n g g + e + + v. for n g < 2. We note that we prefer to use helicity amplitudes for reasons of speed 
and efficiency of the computation. While for top quark decay this is not a critical issue - the number 
of diagrams is sufficiently small even at NNLO - it becomes essential if one considers applicability 
of similar techniques to more complicated processes. Some amplitudes that we present below have 
already appeared in the literature [7] while others are new. 

To deal with the massive top quark, we follow the standard procedure that is based on the 
decomposition of a time- like four-momentum into a sum of two light-like four- momenta [57]. To 
accomplish this, we choose an arbitrary light-like momentum r\ and write 



Pt=Pi-\ r], (A.l) 



D Ptri 



where Sy- = 2pi ■ pj . We note that p\ = by construction. We can now write the top quark spinor as 

U-( Pt ) = (fa + m t )\r,) T L = \l] + \r ] ) * 



M> W (A.2) 

1 777, ^ v 1 

u+(pt) = {h + TO *)Mj7^j = I 1 ) + M t^t, 

where we use the standard notation for spinor-helicity variables. In particular, for pf = 0, \i) = \pi) = 
u+(pi) and \i] = \pi] = U—(pi). It is convenient to choose the auxiliary momentum r\ to be the positron 
momentum p$ to compute helicity amplitudes. 
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We now present results for helicity amplitudes starting from the tree-level process t(p t ) — > b(p-2) + 
v(p§) + e + {p§). We introduce the Breit-Wigner denominator Dw{p\v) = (j>\y ~ m w + imw^w) -1 
and write the amplitude as 

■Ah t (jPt,P2;Ps,Ps) = gwD w {s56)Si 1 i 2 A ht (l,2;5,6) = PS^A^l, 2; 5, 6). (A. 3) 

The amplitude Ah t (1, 2; 5, 6) is the color-stripped amplitude and it only depends on the helicity of the 
top quark h t since helicities of all masslcss fermions are fixed thanks to the left-handed nature of the 
weak current. We find 

A+(l,2;5,6) = 0, A_(l,2;5,6) = (25) [61] . (A.4) 
In this case the sum over colors and helicities is trivial; the amplitude squared reads 

]T A ht ( Pt , P2 ;p 5 ,p 6 ) = \T\ 2 N c \A_(l,2;5,6)\ 2 . (A.5) 

col.hcl 

Next, we consider amplitudes for the process t(pt) — > 6(^2) + gips) + V {P5) + s + (pe). We write it 

as 

■Ah t h g (pt,P2,P3;P5,P6) = SsFT^Ahth^l, 2, 3; 5, 6), (A.6) 

where T a is an SU(3) color matrix normalized as Tr(T a T b ) = 5 ab and h t . g are helicities of the top 
quark and of the gluon, respectively We obtain the following results for the color-stripped helicity 
amplitudes 



A__(l,2,3;5,6) 
A_+(l,2,3;5,6) 



[61](3|p 56 |2](25) [61] (35) 



(3|p t |3][32] [32] ' 

[61](2|p 56 |3](25) [31][63](2|5) 

<3|p t |3](23) (3|fc|3] ' (A . 7) 
(25) [63] 2 



A ++ (l 1 2 ) 3;5,6) = -m t[61](3| ^ |3] , 
A+- = 0, 

where p^e = p§ + pe- The sum over colors is trivial also in this case. We find 

£ \A\ 2 = |-F| 2 <? s 2 (iV c 2 - 1)|A| 2 = 2N c C F \F\ 2 g 2 s \A\ 2 - (A.8) 

col 

The last set of tree amplitudes that we require for our calculation are the ones that describe the 
process with two gluons in the final state t{jpt) ~^ b(p 2 ) + g{pz) + g(Pi) + vips) + s + (p&)- We write 

A hthgs h g4 (pt,P2,P3,P4;p 5 , Pe ) = Tg 2 C 5 ™" 3 ■ Ta ' l )i^ A tHh g<r3 h gai (1, 2, <r(3), <r(4); 5, 6), (A.9) 
where a denotes two permutations of the gluons 53, 54. The color-stripped helicity amplitudes A(l,2,3,4;5,6) 
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are computed to be 

[61] (25) [2\p t p M \2] [61] (4|p 56 |2] (25) (3|p t |2] 



S t34 [32][42][43] s t34 [32][42][3|p t |3) 
[61](45)(3|p t |2] [61] (45) [61] (35) 
[32][42][3|pt|3) [32] [43] [42] [43] ' 

( [61] C_MM! + ( 3 ip*i 4 ] 2 , <; 

" + [ 0> ^43] (34) ^ S234 (2|4) + S t34 (3|p t |3] + (: 
[4|l][6|4](3|p t |4] \ [6|1](2|3) 2 (3|5) 
S t34 [4|3](3|p t |3]y s 23 4(2|4)(3|4) ' 

61](25) /[32](4|p i |3] [32] 2 (24) , (4|p,|3] 5 



[43] (34) V[42](3|p t |3] s 234 [42] s t34 (3|p 4 |3], 
[61](45) Z[32](24) (4|p t |3]\ , [63][31](25) / (4|p t |3] [32] \ , [31][63](45) 



[42] (34) V s 234 (3|pt|3]y [43](3|pt|3] V S tu [42] J [42](3|p t |3] ' 
(25) /[61](2|p 56 |4](4|p t |3] [31][63](2|p 56 |4] , [61](2|p t |3] , 



and 



s t34 V (24)(34)(3| Pt |3] (24)(3|p t |3] (24) (34) 

[41][64](4|p t |3] [31] [43] [64] [31] [63] (23) 2[31][64] 
(34)(3|p t |3] <3|p t |3] + (24)(34) + (34) 



A+— =0, 

[64] 2 (25)(3|p t |4] 

A+-+ = m t - 



s t34 [43][61](3|p t |3]' 

[63] 2 / (25)(4|p t |3] (45) [32](25)\ (A.10) 



[61](3|p t |3] V s t34 [43] + [42] [42] [43] J 



(25) /[43](l|p 34 |6] , [63](2|p 34 |6] 



(34)(3|p t |3] V 5 t34 [61] (24) 

with S( 34 = s 34 — s 3t — s 4t . The color-summed amplitude now reads 



£|^| 2 = 4iV c |.F| 2 5 4 C7 F 



col 



C7 F (|A(l,2,3,4;5,6)| 2 +|^(l,2,4,3;5,6)| i 



(A.11) 



Cf ~ IT I 

B Calculation of one-loop integrals for the soft limit 

In this Appendix, we give an example of how one-loop eikonal integrals needed for the soft current at 
one-loop can be computed. Consider the four-point function, defined in Eq.(4.11) 

d d k 1 mH . 

(B.l) 



i{2n) d (k - p 3 ) 2 fc 2 (2p 6 • k)(2 Pt ■ k - 2p t ■ p 3 ) ' 



where p\ = ntf, p 2 = P 3 = 0. We assume that each propagator contains a term +i0 which defines the 
proper analytic continuation in cases of branch cuts. 
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To compute I4, we first combine k 2 and (k — p 3 ) 2 using Feynman parameters 

(B.2) 



1 r dx 



k 2 (k-p 3 ) 2 J [(k-x P3 ) 2 } 2 


Next, we exponentiate the eikonal quark propagators introducing the integrals over "proper time 

1 



(2 Pt ■k-2pfp 3 + i0)(2fc • Pb + iO) 



J dT 1 dT 2 e t2k - pbTl+l( - 2pt - k - 2pt - p3)T2 . (B.3) 
b 



We now combine Eqs.(B.2, B.3) and shift the loop momentum k — > k + p 3 x. The resulting integral 
becomes 

00 1 

h = i [ dTjdT-2 / dxe li ~ 2pt - p3il - x ^ +2p3 - pbXT ^ [ ^^-l-e^ 2ptT2+2p ^\ (B.4) 

J J J {^r [k z \ z 



The integral over k is easily evaluated 

J (27T)d[k 2 } 2 ~ (4^/2 ^2+^* P^nj . (tf.iJj 

To proceed further, we change variables n = T2^, < ^ < oo and integrate over T2 and x. We obtain 

oo 

r(-e)r(l + 2e) /■ dg (m 2 + 2 Pt -p b ^Y ( e 2 ™ 1 \ 



{Airy/ 2 J {2p t -p 3 + 2p b -p 3 £) \(2 Pb ■ p 3 1+2i (2 P fP 3 ) 1+ 
o 

The final integration over £ is easy to perform in terms of hypergeometric functions. The result is 
shown in Eq.(4.12). 
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